#include "retina.h"
#define PI2 6.28318530717958647692
//***************************************************************************
//******************  The Retina functions   ********************************
#include <CImg.h>

Retina::Retina(double stepp)
{
  step=stepp; input_luminosity_range=255.0; pixels_per_degree=1.0;
  sizeX=0; sizeY=0;
  //Just for safety, really...
  radScheme=0; outerPlexiformLayer=0; contrastGainControl=0; AmacrineCells=0;
  done=false; adjusted=false;
}


void Retina::xmlize()
{
  // General parameters
  add_param(step,"temporal-step__sec");
  add_param( input_luminosity_range , "input-luminosity-range");
  add_param( pixels_per_degree , "pixels-per-degree");

  // Log-polar scheme
  add_child(radScheme,"");
  xmlParam::Register<1, ret_LogPolarScheme>(radScheme , "log-polar-scheme");

  // To be applied right after
  // OPL schemes
  add_child(outerPlexiformLayer,"outer-plexiform-layer");
  xmlParam::Register<1, LinearOPL > ( outerPlexiformLayer, "linear-version" );
  xmlParam::Register_init<1> (outerPlexiformLayer,"linear-version", &Retina::setFather, *this );
  xmlParam::Register<1, GainControlOPL > (outerPlexiformLayer, "luminance-gain-control-version" );
  xmlParam::Register_init<1> (outerPlexiformLayer,"luminance-gain-control-version", &Retina::setFather, *this );

  // Gain control scheme
  add_child ( contrastGainControl, "" );
  xmlParam::Register<1, FastGainControl >( contrastGainControl , "contrast-gain-control" );
  xmlParam::Register_init<1>(contrastGainControl, "contrast-gain-control", &Retina::setFather, *this );
  xmlParam::Register<1, FastGainControl_custom >( contrastGainControl , "contrast-gain-control-custom" );
  xmlParam::Register_init<1>(contrastGainControl, "contrast-gain-control-custom", &Retina::setFather, *this );

  // Directive-Selection Amacrine cells layers
  add_child ( AmacrineCells, "");
  xmlParam::Register<1, AmacrineLayer >( AmacrineCells , "amacrine-layer" );
  xmlParam::Register_init<1>(AmacrineCells, "amacrine-layer", &Retina::setFather, *this );

  // Ganglion layers
  add_child ( ganglionLayers, "" );
  xmlParam::Register<2, GanglionLayer> (ganglionLayers, xmlParam::names ("" , "ganglion-layer") );
  xmlParam::Register_init<2>(ganglionLayers, "ganglion-layer", &Retina::setFather, *this );

  // Directive-Selection Ganglion cells layers
  add_child ( dsganglionLayers, "" );
  xmlParam::Register<2, DSGanglionLayer> (dsganglionLayers, xmlParam::names ("" , "ds-ganglion-layer") );
  xmlParam::Register_init<2>(dsganglionLayers, "ds-ganglion-layer", &Retina::setFather, *this );

}


Retina::~Retina()
{
  delete radScheme;
  delete outerPlexiformLayer;
  delete contrastGainControl;
  delete AmacrineCells;
}


// Return minimum size that an input image should have, given the blur parameters
// required for the retina in its successive layers, for the filtering to be exact
// at least in central pixel of the retina.
int Retina::get_minimumSize()
{
  double sumSigmaSquare=0.0;
  for(uint chann=0;chann<ganglionLayers.size();++chann)
    sumSigmaSquare=max(sumSigmaSquare,
                       pow(ganglionLayers[chann]->sigmaPool,2.0));
  if(outerPlexiformLayer)
    sumSigmaSquare+=(pow(outerPlexiformLayer->sigmaCenter,2.0)
                    +pow(outerPlexiformLayer->sigmaSurround,2.0));
  if(contrastGainControl)
    sumSigmaSquare+=pow(contrastGainControl->sigmaSurround,2.0);

  //taking 3.5*sigma as the 'zero' boundary for Gaussian kernels
  return (int)ceil(7*sqrt(sumSigmaSquare)*pixels_per_degree);
}


void Retina::setFather(RetinaStage *retStage)
{
  retStage->fatherRetina=this;
}


// Helper function:
// Adjust ONE size of ONE spiking channel as a function of input image dimension.
// output code: 0=no prob, 1=auto, 2=too big
int Retina::adjustChannelSize(ret_SpikingChannel* chan, int imSize, const char dim)
{
  if(!chan->sampScheme)
  {
    cerr<<"No proper sampling array ('circular-array' or 'square-array') "
          "has been found in this spiking channel. The retina cannot "
          "initialize its elements."<<endl;
    return(0);
  }
  double *chanSize;
  if(dim=='X') // X
    chanSize=chan->sampScheme->sizeX;
  else // Y
    chanSize=chan->sampScheme->sizeY;

  int ret=0; // "do not touch"
  if( round(*chanSize*pixels_per_degree) > imSize )
    ret=2;   // "need to reduce"
  else if( *chanSize == 0 )
    ret=1;   // "need to set size automatically"

  if( ret )
    *chanSize = (double)imSize/pixels_per_degree; // update channel size
  return ret ;
}


// Adjusting sizeX and sizeY, that is, the size of all maps that are going to be
// allocated by sub-objects of the retina, as a function of the sizes of the input
// images that are going to be fed to the retina...
void Retina::adjustSizes(int imX, int imY)
{
  if(adjusted)
  {
    cerr<<"ERROR in Retina::adjustSizes(): Trying to adjust sizes for the second time!"<<endl;
    exit(0);
  }

  //(1) Determine a 'confidence' zone in the image, without edge effects in the filtering
  int minSizeForBlurs=get_minimumSize();
  if(imX<minSizeForBlurs || imY<minSizeForBlurs)
    cerr<<"WARNING in Retina::adjustSizes(): Given the spatial filtering parameters "
      "(\'xx-sigma-xx\') required for the retina in its successive layers, the input "
      "image should be at least of size "<<minSizeForBlurs<<", to avoid edge effects "
      "in the filtering at central pixel (at least)."<<endl;
  int imX_conf=max(1,imX-minSizeForBlurs+1);
  int imY_conf=max(1,imY-minSizeForBlurs+1);


  //(2) FIX sizeX and sizeY from the spiking channels' size specifications
  int nbSC=0;

  for(uint chann=0;chann<ganglionLayers.size();++chann)
  {
    if(ganglionLayers[chann]->spikingChannel)
    {
      nbSC++;

      // Adjust channel sizes. Because circular arrays have only *one* diameter,
      // it is important to first adjust to the smallest dimension :
      int retX, retY;
      if(imX_conf<imY_conf)
      {
        retX=adjustChannelSize( ganglionLayers[chann]->spikingChannel, imX_conf, 'X');
        retY=adjustChannelSize( ganglionLayers[chann]->spikingChannel, imY_conf, 'Y');
      }
      else
      {
        retY=adjustChannelSize( ganglionLayers[chann]->spikingChannel, imY_conf, 'Y');
        retX=adjustChannelSize( ganglionLayers[chann]->spikingChannel, imX_conf, 'X');
      }

      // Update corresponding retinal sizes
      sizeX=max(sizeX,
               (int)round(*(ganglionLayers[chann]->spikingChannel->sampScheme->sizeX)*pixels_per_degree));
      sizeY=max(sizeY,
               (int)round(*(ganglionLayers[chann]->spikingChannel->sampScheme->sizeY)*pixels_per_degree));

      //Always keep the customer informed.
      int status=max(retX,retY);
      if ( status==1 )
        cout<<"From Retina::adjustSizes(): A spiking array <"
        << xmlParam::last_xmlName( *(ganglionLayers[chann]->spikingChannel->sampScheme) ) // careful using this function :p
        <<"> is being automatically set to fit the size of input image (minus a border for edge effects)."<<endl;
      else if ( status==2 )
        cout<<"WARNING from Retina::adjustSizes(): A spiking array <"
        << xmlParam::last_xmlName( *(ganglionLayers[chann]->spikingChannel->sampScheme) ) // careful using this function :p
        <<"> is being automatically REDUCED to fit the size of input image (minus a border for edge effects)."<<endl;
    }
  }

  //(2bis) If the retina has no spiking channel, then sizeX and sizeY are just defined by the input image!
  if(!nbSC)
    sizeX=imX_conf, sizeY=imY_conf;

  //(3) add back the 'border for blur' to the retinal maps, keeping the constraint that sizeX<=imX
  sizeX=min(sizeX+minSizeForBlurs-1, imX);
  sizeY=min(sizeY+minSizeForBlurs-1, imY);
  adjusted=true;
}



void Retina::allocateValues(bool verbose, double initValue)
{
  if(!adjusted)
  {
    cerr<<"ERROR in Retina::allocateValues(): use adjustSizes(imX, imY) first!"<<endl;
    exit(0);
  }

  double initVal = (initValue==HUGE_VAL)?
                   input_luminosity_range / 2.0 :
                   initValue ;

  if(verbose)
  {
    cout<<endl<<"The retina is currently allocating all its elements!"<<endl;
    cout<<"Size of its maps was fixed at: ("<<sizeX<<","<<sizeY<<")."<<endl;
    cout<<"Initialization of the retina (OPL only) by a uniform screen"
      " of luminance "<<initVal<<"."<<endl;

    cout<<endl<<"(1) OUTER PLEXIFORM LAYER:"<<endl;
  }

  if(outerPlexiformLayer)
    outerPlexiformLayer->allocateValues(verbose, initVal );
  else
  {
    cerr<<"ERROR: Your retina does not have an Outer Plexiform Layer, "
      "which is its minimal element!"<<endl;
    exit(0);
  }

  if(verbose)
    cout<<endl<<"(2) CONTRAST GAIN CONTROL MECHANISM:"<<endl;

  if(contrastGainControl)
    contrastGainControl->allocateValues(verbose);
  else if(verbose)
      cout<<"No Contrast Gain Control mechanism was defined."<<endl;

  if(verbose)
    cout<<endl<<"(3) AMACRINE CELLS LAYER:"<<endl;

  if(AmacrineCells)
      AmacrineCells->allocateValues(verbose);

  else if(verbose)
      cout<<"No Amacrine Cells Layer was defined."<<endl;

  if(verbose)
    cout<<endl<<"(4) SPIKING LAYERS OF GANGLION CELLS:"<<endl;

  if(ganglionLayers.size())
  {
    int nbUnits=0, nbSpkCh=0, nbCh=ganglionLayers.size();
    for(int chann=0;chann<nbCh;++chann)
    {
      ganglionLayers[chann]->allocateValues(verbose);
      if(ganglionLayers[chann]->spikingChannel)
      {
        nbSpkCh++;
        nbUnits+=ganglionLayers[chann]->spikingChannel->units.size();
      }
    }
    if(verbose)
      cout<<nbCh<<" layers of ganglion cells were defined, amongst which "
        <<nbSpkCh<<" are spiking channels, for a total of "<<nbUnits
        <<" spiking ganglion cells."<<endl;
  }
  else if(verbose)
      cout<<"No layer of ganglion cells were defined. Retinal output will be "
        "an analog spatio-temporal map corresponding to the potential of a "
        "layer of bipolar cells."<<endl;

  if(verbose)
    cout<<endl<<"(5) SPIKING LAYERS OF DIRECTION-SELECTIVE GANGLION CELLS:"<<endl;

  if(dsganglionLayers.size())
  {
      for(int chann=0; chann<dsganglionLayers.size(); ++chann){
          dsganglionLayers[chann]->allocateValues(verbose,chann);
      }
      if(verbose)
      cout<<"Size of the DSGC array:" <<dsganglionLayers.size()<<endl;
  }

  else if(verbose)
      cout<<"No DSGC Layer was defined."<<endl;

  done=true;
}

void Retina::feedInput(const CImg<double>& image)
{
  if(!done)
  {
    adjustSizes(image.width(),image.height());
    allocateValues();
  }
  //if the input image should be bigger than necessary, reduce it:
  if((image.width()==sizeX)&(image.height()==sizeY))
  {
    input.const_targets=&image;
    input.targets=0;
  }
  else
  {
    stockInput = image.get_crop((image.width()-sizeX)/2,
                                (image.height()-sizeY)/2,
                                (image.width()-sizeX)/2+sizeX-1,
                                (image.height()-sizeY)/2+sizeY-1);
    input.targets=&stockInput;
  }
  outerPlexiformLayer->feedInput(input.read());
  if(contrastGainControl)
    contrastGainControl->feedInput(outerPlexiformLayer->outputSignal->read());
  if(AmacrineCells)
      if(contrastGainControl)
          AmacrineCells->feedInput(contrastGainControl->bipolarSignal->read());
      else
          AmacrineCells->feedInput(outerPlexiformLayer->outputSignal->read());
  for(uint chann=0;chann<ganglionLayers.size();++chann)
    if(contrastGainControl)
      ganglionLayers[chann]->feedInput(contrastGainControl->bipolarSignal->read());
    else
      ganglionLayers[chann]->feedInput(outerPlexiformLayer->outputSignal->read());
  for(uint chann=0; chann<dsganglionLayers.size();++chann)
      if(contrastGainControl)
          dsganglionLayers[chann]->feedInput(contrastGainControl->bipolarSignal->read());
      else
          dsganglionLayers[chann]->feedInput(outerPlexiformLayer->outputSignal->read());
}

void Retina::tempStep()
{
  if(!done)
    cerr<<"ERROR in Retina::tempStep(): feedInput first!"<<endl;
  outerPlexiformLayer->tempStep();
  if(contrastGainControl)
    contrastGainControl->tempStep();
  if(AmacrineCells)
      AmacrineCells->tempStep();
  for(uint chann=0;chann<ganglionLayers.size();++chann)
    ganglionLayers[chann]->tempStep();
  for(uint chann=0;chann<dsganglionLayers.size();++chann)
    dsganglionLayers[chann]->tempStep();
}



//***************************************************************************
//****************** Base class RetinaStage *********************************

void RetinaStage::checkFather()
{
  if( !fatherRetina )
  {
    cerr<<"ERROR in this retinal stage: Its fatherRetina has not been set."<<endl;
    exit(0);
  }
}


//****************************************************************************
//******** Functions for OPL classes LinearOPL and GainControlOPL   **********



//******************LinearOPL ******************:

LinearOPL::LinearOPL()
{
  sigmaCenter=0.0;
  tauCenter=0.01;
  nCenter=0;
  sigmaSurround=1.0;
  tauSurround=0.01;
  wSurround=0.0;
  leakyHeatEq=0;
  ampOPL=1.0; relative_ampOPL=1.0;

  excitCells = new RetinaMapFilter();
  inhibCells = new RetinaMapFilter();
  undershoot = 0 ;
  outputSignal=0, receptorSignal=0, horizontalSignal=0;//, toAmacrineSignal=0;
  //Vbip.fill(0);

}

void LinearOPL::xmlize()
{
  add_param(sigmaCenter,"center-sigma__deg");
  add_param(tauCenter,"center-tau__sec");
  add_param(nCenter,"center-n__uint");
  add_param(sigmaSurround,"surround-sigma__deg");
  add_param(tauSurround,"surround-tau__sec");
  add_param(wSurround,"opl-relative-weight");
  add_param(leakyHeatEq,"leaky-heat-equation");
  add_param(relative_ampOPL,"opl-amplification");

  //The undershoot is now added as a subobject:
  add_child( undershoot , "" ) ;
  xmlParam::Register<1,LinearOPL::Undershoot> ( undershoot , "undershoot" ) ;
}

void LinearOPL::Undershoot::xmlize()
{
  add_param(wAdap,"relative-weight");
  add_param(tauAdap,"tau__sec");
}


LinearOPL::~LinearOPL()
{
  delete excitCells;
  delete inhibCells;
  delete outputSignal;
  delete undershoot ;
}


void LinearOPL::allocateValues( bool verbose, double initValue)
{
  if(verbose)
  {
    cout<<"Allocating values for the LINEAR version LinearOPL."<<endl;
    if(leakyHeatEq)
      cout<<"Low-pass filter for each layer of cells is a non-separable "
        "leaky heat equation."<<endl;
    else
      cout<<"Low-pass filter for each layer of cells is a separable "
        "spatial_gaussian / temporal_negative_exponential."<<endl;
    if(undershoot)
      cout<<"A linear undershoot (transient) has been added."<<endl;
  }

  checkFather();
  int sx=fatherRetina->sizeX;
  int sy=fatherRetina->sizeY;
  double st=fatherRetina->step;
  double ppd=fatherRetina->pixels_per_degree;
  ret_LogPolarScheme *rs=fatherRetina->radScheme;

  //relative_ampOPL must have values of order [1-10] for the OPL output to have typical values of order 1.
  if ( (relative_ampOPL<1) || (relative_ampOPL>20) )
    cerr<<"WARNING from the LinearOPL object: The parameter LinearOPL::opl-amplification "
      "does not seem well suited. It is a RELATIVE factor with respect to input "
      "luminosity on the retina. It should always be set between 1 and 20 for the "
      "output range of the OPL to have values of order 1. CHANGE IT in your retina "
      "definition file!"<<endl;
  ampOPL = relative_ampOPL / fatherRetina->input_luminosity_range ;

  excitCells->set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
  if(leakyHeatEq)
  {
    excitCells -> leakyHeatEquation( tauCenter , sigmaCenter );
    if ( nCenter != 0 )
      cerr<<"WARNING from the LinearOPL object: You have chosen a 'leaky heat equation' "
        "filter, so parameter center-n__uint is inefficient ! It will be considered "
        "as having value zero (simple exponential filter)."<<endl;
  }
  else
  {
    excitCells -> ExpCascade(tauCenter , nCenter);
    excitCells -> set_sigmaPool(sigmaCenter);
  }
  (*excitCells) *= ampOPL;

  //We set the initial input of the filters at half the retinal input_luminosity_range
  excitCells -> set_initial_input_value ( initValue );
  excitCells -> allocateValues();
  receptorSignal = excitCells;


  if(undershoot)
  {
    BaseRecFilter adapTr( st );
    adapTr.Exp ( undershoot->tauAdap ).warningVerbose=false;
    adapTr*= undershoot->wAdap ;

    (*undershoot) .set_sizeX(sx).set_sizeY(sy).set_step(st) ;
    (*undershoot) -= adapTr ;

    undershoot -> set_initial_input_value ( excitCells->read()(sx/2,sy/2) );
    undershoot -> allocateValues();
    receptorSignal = undershoot ;
  }


  if(wSurround)
  {
    inhibCells -> set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
    if(leakyHeatEq)
      inhibCells -> leakyHeatEquation( tauSurround , sigmaSurround );
    else
    {
      inhibCells -> Exp(tauSurround);
      inhibCells -> set_sigmaPool(sigmaSurround);
    }
    (*inhibCells) *= wSurround;

    //We set the initial input of the filters at half the retinal input_luminosity_range
    //NOTA : receptorSignal = undershoot if it is present, and = excitCells otherwise.
    inhibCells -> set_initial_input_value ( receptorSignal->read()(sx/2,sy/2) );
    inhibCells -> allocateValues();
    horizontalSignal = inhibCells;

    stock_output = CImg<double>(sx,sy,1,1,0.0);
    outputSignal = new CImgReader;
    outputSignal -> targets= &stock_output;
  }
  else
    outputSignal = receptorSignal ;

}


void LinearOPL::feedInput(const CImg<double>& image)
{
  excitCells -> feedInput( image );
  if( undershoot )
    undershoot -> feedInput ( excitCells -> read() );
  //NOTA : receptorSignal = undershoot if it is present, and = excitCells otherwise.
  if( wSurround )
    inhibCells -> feedInput( receptorSignal -> read() );
}

void LinearOPL::tempStep()
{
  excitCells -> tempStep();
  if ( undershoot )
    undershoot -> tempStep();
  if(wSurround)
  {
    inhibCells -> tempStep();
    stock_output = receptorSignal -> read() - inhibCells -> read();
  }

}



//****************** GainControlOPL ******************:

GainControlOPL::GainControlOPL()
{
  ampPhoto=1.0;
  sigmaCenter=0.0;
  tauPhoto=0.01;
  sigmaSurround=1.0;
  tauSurround=0.01;
  wSurround=0.0;
  gLeak=10.0;
  ampFeedback=200.0;
  nernstFeedback=0.0;

  phototransFilter=RetinaMapFilter();
  excitCells=CondCellMap();
  controlCond=RetinaMapFilter();
  outputSignal=0, receptorSignal=0, horizontalSignal=0;
}

void GainControlOPL::xmlize()
{
  add_param(ampPhoto,"input-amplification-phototransduction");
  add_param(sigmaCenter,"sigma-receptors__deg");
  add_param(tauPhoto,"tau-phototransduction__sec");
  add_param(sigmaSurround,"sigma-horizontal-cells__deg");
  add_param(tauSurround,"additional-tau-horizontal-cells__sec");
  add_param(wSurround,"relative-weight-surround-center");
  add_param(gLeak,"inert-leak-in-receptors");
  add_param(ampFeedback,"horizontal-feedback-amplification");
  add_param(nernstFeedback,"horizontal-feedback-nernst-potential");
}


GainControlOPL::~GainControlOPL()
{
  if((wSurround)&&(outputSignal)) delete outputSignal;
}



void GainControlOPL::allocateValues( bool verbose, double initValue)
{
  if(verbose)
  {
    cout<<"Allocating values for a CONDUCTANCE-BASED OPL (luminance gain-control)."<<endl;
    cout<<"Very experimental element. the initValue method does not work (amongst other things!)"<<endl;
  }

  checkFather();
  int sx=fatherRetina->sizeX;
  int sy=fatherRetina->sizeY;
  double st=fatherRetina->step;
  double ppd=fatherRetina->pixels_per_degree;
  ret_LogPolarScheme *rs=fatherRetina->radScheme;

  phototransFilter.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
  phototransFilter.Exp(tauPhoto);
  phototransFilter.set_sigmaPool(sigmaCenter);
  phototransFilter*=ampPhoto;

  excitCells.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
  excitCells.set_g_leak(gLeak);
  int index=excitCells.new_ionicChannel(nernstFeedback);
  excitCells.inputSynapse(index).linearAmplification(ampFeedback);
  excitCells.allocateValues();
  receptorSignal=&excitCells;

  controlCond.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
  controlCond.Exp(tauSurround); controlCond.set_sigmaPool(sigmaSurround);
  controlCond.inputSynapse().rectification(0.0);
  controlCond.allocateValues();
  horizontalSignal=&controlCond;

  if(wSurround)
  {
    stock_output=CImg<double>(sx,sy);
    outputSignal=new CImgReader;
    outputSignal->targets=&stock_output;
  }
  else
    outputSignal=&excitCells;


}

void GainControlOPL::feedInput(const CImg<double>& image)
{
  excitCells.feedCurrent(image);//the input current corresponds to port 0...
  excitCells.feedCond(controlCond.read(),1,true);//the feedback conductance corresponds to port 1...
  controlCond.feedInput(excitCells.read());
}

void GainControlOPL::tempStep()
{
  excitCells.tempStep();
  controlCond.tempStep();
  if(wSurround)
    stock_output=excitCells.read()-wSurround*controlCond.read();
}




// ******************************************************
// ****************** FastGainControl:  *****************
// ******************************************************



FastGainControl::FastGainControl() :  	vThresh(0.0), Einh(0.0), exponent(2.0), symmetrize(1),
    inhibCells(0), wTransient(0), tauTransient(0.03)
{
  ampInputCurrent=100.0;
  sigmaSurround=1.0;
  tauSurround=0.01;
  gLeak=10.0;
  ampFeedback=100.0;

  excitCells=CondCellMap();
  controlCond=RetinaMapFilter();
  bipolarSignal=0, adaptationSignal=0;
}


void FastGainControl::xmlize()
{
  add_param(ampInputCurrent,"opl-amplification__Hz");
  add_param(gLeak,"bipolar-inert-leaks__Hz");

  add_param(sigmaSurround,"adaptation-sigma__deg");
  add_param(tauSurround,"adaptation-tau__sec");
  add_param(ampFeedback,"adaptation-feedback-amplification__Hz");
}


void FastGainControl::allocateValues( bool verbose)
{

  if(verbose)
    cout<<"Allocating values for a conductance-based contrast-gain-control."<<endl;

  checkFather();
  int sx=fatherRetina->sizeX;
  int sy=fatherRetina->sizeY;
  double st=fatherRetina->step;
  double ppd=fatherRetina->pixels_per_degree;
  ret_LogPolarScheme *rs=fatherRetina->radScheme;

  // *********** (1) LAYER OF BIPOLAR CELLS:

  excitCells.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
  excitCells.set_g_leak(gLeak);

  // (a) Input current in port number 0.
  excitCells.inputSynapse().linearAmplification(ampInputCurrent);

  // (b) Adapting inhibitory (or leak, if Einh=0) conductance in port number 1 -- providing CGC.
  // Einh=0 should be the law if an ON/OFF symmetry is present.
  excitCells.new_ionicChannel( Einh );

  // (c) Optional linear transient, if wTransient!=0 :
  if(wTransient)
    excitCells.new_currentInput();

  // (d) Allocating the retinal map. Pointing the reader on the map.
  excitCells.allocateValues();
  bipolarSignal=&excitCells;



  // *************** (2) ADAPTING CONDUCTANCES IN THE MEMBRANE

  controlCond.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);

  // (a) Synaptic transmission ==> Parameters for this stage can only be fixed through the 'CUSTOM' version.

  controlCond.inputSynapse().rectification( vThresh , ampFeedback );
  controlCond.inputSynapse().set_power_exp ( exponent ) ;
  if (symmetrize)
  {
    controlCond.new_currentInput();
    controlCond.inputSynapse(1).rectification( -vThresh , -ampFeedback ); //The simultaneous 'OFF'-bipolar-control
    controlCond.inputSynapse(1).set_power_exp ( exponent ) ;
  }

  // (b) Spatio-temporal averaging properties

  controlCond.Exp(tauSurround);
  controlCond.set_sigmaPool( sigmaSurround );

  // (c) Allocating the retinal map. Pointing the reader on the map.
  controlCond.allocateValues();
  adaptationSignal=&controlCond;


  // ***************** (3) OPTIONAL ADDTIONAL LINEAR HIGH-PASS
  if(wTransient)
  {
    inhibCells = new RetinaMapFilter() ;
    inhibCells -> set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
    inhibCells -> Exp(tauTransient) ;
    (*inhibCells) *= wTransient ;
    inhibCells -> allocateValues() ;
  }

}

void FastGainControl::feedInput(const CImg<double>& image)
{
  // *******  (1) Inputs to bipolar cells

  // input current in port 0:
  excitCells.feedCurrent(image);
  //feedback conductance corresponds to port 1:
  excitCells.feedCond(controlCond.read(),1,true);
  //Optionally, linear feedback corresponds to port 2:
  if(wTransient)
    excitCells.feedCurrent(inhibCells->read(),2);

  // *******  (2) Inputs to adapting conductance:
  controlCond.feedInput(excitCells.read(),0);
  if(symmetrize)
    controlCond.feedInput(excitCells.read(),1); //Input from 'ON' and 'OFF' bipolar cells.

  // *******  (3) OPTIONALLY, Input to the supplementary transient (short-range amacrine cells)
  if (inhibCells)
    inhibCells->feedInput( excitCells.read() );
}

void FastGainControl::tempStep()
{
  excitCells.tempStep() ;
  controlCond.tempStep() ;
  if(inhibCells)
    inhibCells->tempStep() ;
}


//****************** FastGainControl_custom ******************:


void FastGainControl_custom::xmlize()
{
  FastGainControl::xmlize();
  add_param(exponent,"adaptation-powerlaw-exponent");
  add_param(Einh,"adaptation-E-inh");
  add_param(symmetrize,"adaptation-ONOFF-symmetry");
  add_param(vThresh,"bipolar-to-adaptation-threshold");

  add_param( wTransient , "transient-relative-weight");
  add_param( tauTransient , "transient-tau__sec");
}


// ********************************************************
// ****************  GanglionLayer  ***********************
// ********************************************************

GanglionLayer::GanglionLayer()
{
  sign=1.0;
  vThresh=0.0;
  valThresh=70.0;
  ampInputCurrent=100.0;
  tauTransient=0.04;
  wTransient=0.75;
  sigmaPool=0.0;

  spikingChannel=0;
}


void GanglionLayer::xmlize()
{
  add_param(sign,"sign");
  add_param(tauTransient,"transient-tau__sec");
  add_param(wTransient,"transient-relative-weight");
  add_param(vThresh,"bipolar-linear-threshold");
  add_param(valThresh,"value-at-linear-threshold__Hz");
  add_param(ampInputCurrent,"bipolar-amplification__Hz");
  add_param(sigmaPool,"sigma-pool__deg");

  // auto-registration of the pointer:
  add_child(spikingChannel,"");
  xmlParam::Register<1,ret_SpikingChannel>(spikingChannel,"spiking-channel");
  xmlParam::Register<1,ret_SpikingChannel_adap>(spikingChannel,"spiking-channel-adap");
}


// xmlParam initialization function (right after full XML load)
void GanglionLayer::xmlinit_after_children()
{
  if(spikingChannel)
  {
    spikingChannel->radScheme=fatherRetina->radScheme;
    if(spikingChannel->sampScheme)
      spikingChannel->sampScheme->father=spikingChannel;
  }
}


GanglionLayer::~GanglionLayer()
{
  delete spikingChannel;
}

void GanglionLayer::allocateValues( bool verbose)
{
  sign=(sign<0)? -1.0:1.0;

  checkFather();
  int sx=fatherRetina->sizeX;
  int sy=fatherRetina->sizeY;
  double st=fatherRetina->step;
  double ppd=fatherRetina->pixels_per_degree;
  ret_LogPolarScheme *rs=fatherRetina->radScheme;

  inputCurrent.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
  supplementaryTransient.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);

  BaseRecFilter amacrineCuttingInGanglionar(st);
  amacrineCuttingInGanglionar.Exp(tauTransient).warningVerbose=false;
  amacrineCuttingInGanglionar*=wTransient;

  //This is what we had before!!: transient AFTER the synaptic rectif/pooling
  //if(wTransient) inputCurrent-=amacrineCuttingInGanglionar;

  //This is what we propose now: transient BEFORE the synaptic rectif/pooling
  if(wTransient)
    supplementaryTransient-=amacrineCuttingInGanglionar;
  supplementaryTransient.allocateValues();

  //Check the corresponding change in functions feedInput() and tempStep() !

  inputCurrent.inputSynapse()
              .rectification(sign*vThresh,sign*ampInputCurrent,valThresh);
  inputCurrent.set_sigmaPool(sigmaPool);
  inputCurrent.allocateValues();

  if(spikingChannel)
  {
    spikingChannel->set_step(st).set_ratiopix(ppd).set_radScheme(rs);
    spikingChannel->allocateValues();
  }
}

void GanglionLayer::feedInput(const CImg<double>& image)
{
  //This is what we had before!!: transient AFTER the synaptic rectif/pooling
  //inputCurrent.feedInput(image);

  //This is what we propose now: transient BEFORE the synaptic rectif/pooling
  if(wTransient)
  {
    supplementaryTransient.feedInput( image );
    inputCurrent.feedInput( supplementaryTransient.read() );
  }
  else
    inputCurrent.feedInput( image );
  if(spikingChannel)
    spikingChannel->feedInput(inputCurrent.read());
}

void GanglionLayer::tempStep()
{
  if(wTransient)
    supplementaryTransient.tempStep();
  inputCurrent.tempStep();
}

// ********************************************************
//******************* Amacrine Layer **********************
// ********************************************************

AmacrineLayer::AmacrineLayer()
{
    sign = 1.0;
    vThresh = 0.0;
    valThresh = 70.0;
    ampInputCurrent = 100.0;
    ampInputCurrentAm = 100.0;
    tauTransient = 0.04;
    wTransient = 0.75;
    sigmaPool = 0.0;

    nbCells = 0;
    density = 0.0;
    directions = 0;
    radius = 10;
    out_proportion = 0.33;
    radiusOutputStarts = radius * (1.0-out_proportion);
    sx = 0;
    sy = 0;
    alpha = 80;
    beta = 6;
    k1 = 0.002;
    k2 = 0.02;
    theta1 = -0.05;
    theta2 = 0.3;
    gclRest = 100.0/72.0;
    gclBound = 100.0/2.4;
    Ecl = -0.06647818;
    outputSignal = 0;
    delta=1.0/3.0;
}

void AmacrineLayer::xmlize()
{
  add_param(tauTransient,"transient-tau__sec");
  add_param(wTransient,"transient-relative-weight");
  add_param(vThresh,"bipolar-linear-threshold");
  add_param(valThresh,"value-at-linear-threshold__Hz");
  add_param(ampInputCurrent,"bipolar-amplification__Hz");
  add_param(ampInputCurrentAm,"amacrine-amplification__Hz");
  add_param(sigmaPool,"sigma-pool__deg");

  add_param(gammaAm,"gamma-dendrite-filter");
  add_param(amLeak,"amacrine-leak");
  add_param(directions,"number-of-directions");
  add_param(nbCells,"number-of-cells");
  add_param(density,"density");
  add_param(radius, "radius");
  add_param(out_proportion, "output-proportion");
  add_param(alpha, "alfa-gcl-conductance");
  add_param(beta, "beta-gcl-conductance");
  add_param(k1, "k1-gcl-conductance");
  add_param(k2, "k2-gcl-conductance");
  add_param(theta1, "theta1-gcl-conductance");
  add_param(theta2, "theta2-gcl-conductance");
  add_param(gclBound, "gcl-bound-conductance");
  add_param(gclRest, "gcl-rest-conductance");
  add_param(Ecl, "clorure-equilibrium-potential");
  add_param(delta,"soma-interaction-ratio");
}

AmacrineLayer::~AmacrineLayer()
{
//    for(int n=0; n<AmCellLayer.size() ; ++n)
//        AmCellLayer.pop_back();
    if(AmCellLayer.size())
        AmCellLayer.clear();
}

void AmacrineLayer::allocateValues(bool verbose)
{
    checkFather();
    if(verbose)
        cout<<"Allocating values for the amacrine cells layer."<<endl;
    sx=fatherRetina->sizeX;
    sy=fatherRetina->sizeY;
    double st=fatherRetina->step;
    double ppd=fatherRetina->pixels_per_degree;
    ret_LogPolarScheme *rs=fatherRetina->radScheme;
    radiusOutputStarts = radius * (1.0-out_proportion);
    double sxx=sx, syy=sy;
    if(density == 0.0)
    density += sqrt(nbCells/(sx * sy));

    if(fatherRetina->radScheme) cout<<"Warning: the Amacrine cells layer and the hole Direction-Selective system has not been made yet for a radial scheme. All the cells and maps are implemented on square maps."<<endl;

    if(((int)directions%2)!=0) cout<<"Warning: the number of directions in the amacrine cells layer isn't an even number, the directions wont be simetrical."<<endl;
    inputCurrent.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
    supplementaryTransient.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);

    BaseRecFilter amacrineCuttingInGanglionar(st);
    amacrineCuttingInGanglionar.Exp(tauTransient).warningVerbose=false;
    amacrineCuttingInGanglionar*=wTransient;
    if(wTransient)
      supplementaryTransient-=amacrineCuttingInGanglionar;
    supplementaryTransient.allocateValues();

    inputCurrent.inputSynapse()
                .rectification(sign*vThresh,sign*ampInputCurrent,valThresh);
    inputCurrent.set_sigmaPool(sigmaPool);
    inputCurrent.allocateValues();

    //Creation of the dendrites as independent neurons and somas (that's why "<=directions", last one are the somas)
    for(int n=0; n<=directions; ++n)
    {
        AmCellLayer.push_back(new AmacrineLayers);
        AmCellLayer.back()->spikingChannel = new ret_SpikingChannel;
        AmCellLayer.back()->spikingChannel->sampScheme = new sqSampling;
        AmCellLayer.back()->spikingChannel->sampScheme->sizeX=&sxx;
        AmCellLayer.back()->spikingChannel->sampScheme->sizeY=&syy;
        AmCellLayer.back()->spikingChannel->sampScheme->Density=&density;
        AmCellLayer.back()->spikingChannel->sampScheme->set();
        AmCellLayer.back()->spikingChannel->sampScheme->father= AmCellLayer.back()->spikingChannel;

        AmCellLayer.back()->spikingChannel->set_step(st).set_ratiopix(ppd).set_radScheme(rs);
        AmCellLayer.back()->spikingChannel->allocateValues();
        for(int m=0; m<AmCellLayer.back()->spikingChannel->units.size();++m)
        {
            AmCellLayer.back()->spikingChannel->units[m]->thres=99999;            // This sets the threshold potential too high. This maintains the normal states of the neurons, which doesn't spike.
            AmCellLayer.back()->spikingChannel->units[m]->set_g_leak(amLeak);
            AmCellLayer.back()->amCondPort = AmCellLayer.back()->spikingChannel->units[m]->new_ionicChannel(Ecl);
            AmCellLayer.back()->spikingChannel->units[m]->sc1 = 0.0;
            AmCellLayer.back()->spikingChannel->units[m]->sc2 = 0.0;
        }
        AmCellLayer.back()->distalComp.resize(AmCellLayer.back()->spikingChannel->units.size());
        AmCellLayer.back()->stock_output = CImg<double>(sx,sy,1,1,0.0);
        AmCellLayer.back()->outputSignal = new CImgReader;
        AmCellLayer.back()->outputSignal -> targets = &(AmCellLayer.back()->stock_output);
        AmCellLayer.back()->stock_input = CImg<double>(sx,sy,1,1,0.0);
        AmCellLayer.back()->inputSignal = new CImgReader;
        AmCellLayer.back()->inputSignal -> targets = &(AmCellLayer.back()->stock_input);
    }
    //Making soma connections with its own dendrites
    for(int n=0; n<AmCellLayer[directions]->spikingChannel->units.size(); n++)
    {
        for(int m=0; m<directions; m++)
            AmCellLayer[directions]->addDistalComp(n, m, n);
    }


    double xPos=0.0, yPos=0.0, xImg=0.0, yImg=0.0, pixel_distance=0.0, pixel_angle=0.0;
    double init_ang=0.0, fin_ang=0.0, pixel_dist=0.0, pixel_ang=0.0;
    double initial_ang = 0.0, final_ang = 0.0, size_ang = 2.0*(M_PI)/directions;
    double direction_ang = 0.0;
    double r2 = radius*radius;
    CImg<double> temp(sx,sy,1,1,0.0);
    bool connected = false;

    //for each cell is necessary to know how many distal compartments has connected
    for(int n=0; n<directions; ++n)
    {
        direction_ang = n*size_ang;

        if(n==0) initial_ang = PI2 - size_ang/2.0;
        else initial_ang = direction_ang - size_ang/2.0;

        final_ang = direction_ang + size_ang/2.0;
        for(int m=0; m<AmCellLayer.back()->spikingChannel->units.size();++m)
        {
            //First of all, must to be related the connection with the soma
            AmCellLayer[n]->addDistalComp(m, directions, m);

            //Then, creates an image with the output section
            temp.fill(0.0);
            xPos = AmCellLayer[n]->spikingChannel->units[m]->xOffset + ( sx/2.0 );
            yPos = ( sy/2.0 ) - AmCellLayer[n]->spikingChannel->units[m]->yOffset;
            for(double x = floor(xPos - radius); x<= xPos+radius; x++)
            {
                for(double y = floor(yPos - radius); y<= yPos+radius; y++)
                {
                    // if the section exceds the image size, continue
                    if(x < 0 || x >= temp.width() || y < 0 || y>= temp.height()) continue;
                    xImg = x - xPos;
                    yImg = - ( y - yPos );
                    pixel_distance = sqrt (pow(xImg,2)+pow(yImg,2));
                    pixel_angle = atan2 (yImg,xImg);
                    if (pixel_angle<0) pixel_angle = PI2 + pixel_angle;

                    if( (pixel_distance>= radiusOutputStarts && pixel_distance<radius)
                        &&  (   ( (n==0) && ( (pixel_angle>=initial_ang) || (pixel_angle<final_ang) ) )
                               || ( (pixel_angle>=initial_ang) && (pixel_angle<final_ang) ) ) )
                    {
                        temp(x,y) = 1.0;
                    }
                }
            }

            //With this image, reviews if this section is connected to another amacrine cell
            for(int j=0; j<directions; j++)
                for(int k=0; k<AmCellLayer.back()->spikingChannel->units.size();++k)
                {
                    connected = false;

                    // Same cell can not be connected
                    // or if the distance^2 is grater than the minimum to be connected (radius^2)
                    if(k == m || ( pow(AmCellLayer[n]->spikingChannel->units[m]->xOffset
                                       -AmCellLayer[j]->spikingChannel->units[k]->xOffset,2) +
                                   pow(AmCellLayer[n]->spikingChannel->units[m]->yOffset
                                       -AmCellLayer[j]->spikingChannel->units[k]->yOffset,2) ) > r2)
                        continue;

                    xPos = AmCellLayer[j]->spikingChannel->units[k]->xOffset + ( sx/2.0 );
                    yPos = ( sy/2.0 ) - AmCellLayer[j]->spikingChannel->units[k]->yOffset;
                    if(j==0) init_ang = PI2 - size_ang/2.0;
                    else init_ang = j*size_ang - size_ang/2.0;
                    fin_ang = j*size_ang + size_ang/2.0;
                    for(double x = xPos - radius; x<= xPos+radius; x++)
                    {
                        for(double y = yPos - radius; y<= yPos+radius; y++)
                        {
                            if(x < 0 || x >= temp.width() || y < 0 || y>= temp.height()) continue;
                            xImg = x - xPos;
                            yImg = - ( y - yPos );
                            pixel_dist = sqrt (pow(xImg,2)+pow(yImg,2));
                            pixel_ang = atan2 (yImg,xImg);
                            if (pixel_ang<0) pixel_ang = PI2 + pixel_ang;
                            if( temp(x,y)==1.0 && (pixel_dist< radius)
                                &&  (   ( (j==0) && ( (pixel_ang>=init_ang) || (pixel_ang<fin_ang) ) )
                                       || ( (pixel_ang>=init_ang) && (pixel_ang<fin_ang) ) ) )
                            {
                                // If is not done yet,
                                // to cell k in layer j is added a connection to cell m in layer n
                                if(!connected)
                                {
                                    AmCellLayer[j]->addDistalComp(k, n, m);
                                    connected = true;
                                }
                                // In case that the cells are already connected, increase the area of the
                                // overlap
                                else
                                    AmCellLayer[j]->distalComp[k].back()->increaseArea();
                            }
                        }
                        //if (cont) break;
                    }
                }
        }
    }

    stock_output = CImg<double>(sx,sy,1,1,0.0);
    outputSignal = new CImgReader;
    outputSignal -> targets = &stock_output;
    nbCells = AmCellLayer[0]->spikingChannel->units.size();
    if(verbose)
    cout<<"Created amacrine cells layer with "<<nbCells<<" cells, "<<directions<<" dendrites directions and a radius of "
       <<radius<<" pxls ("<<radius/ppd<<" retinal degrees)."<<endl;

    char outfile[200];
    snprintf(outfile,200,"../../../workspace/VR/VR_data_c%d_d%d_r%d.txt", (int)nbCells,(int)directions,(int)radius);
    cout<<outfile<<endl;
    file.open(outfile);
    for(int i=0; i<directions; ++i )
    {
        snprintf(outfile,200,"../../../workspace/VR/SC2_%d.txt", i);
        //cout<<outfile<<endl;
        ofstream f;
        f.open(outfile);
        f.close();
    }
    for(int i=0; i<directions; ++i )
    {
        snprintf(outfile,200,"../../../workspace/VR/SC1_%d.txt", i);
        //cout<<outfile<<endl;
        ofstream f;
        f.open(outfile);
        f.close();
    }

    /*for(int i=1; i<=nbCells; i++)
        for(int j=1; j<=nbCells; j++)
        file<<"Sm"<<i<<"_in_C"<<j<<"\t";
    for(int i=1; i<=directions; i++)
        for(int j=1; j<=nbCells; j++)
            file<<"D"<<i<<"C"<<j<<"_Ibip\tD"<<i<<"C"<<j<<"_Iam\tD"<<i<<"C"<<j<<"_Is\t";
    for(int i=1; i<=directions; i++)
        for(int j=1; j<=nbCells; j++)
            file<<"Vd"<<i<<"c"<<j<<"\t";
    file<<endl;


    file<<"somaV\t";
    for(int j=1; j<=directions; j++)
        file<<"D"<<j<<"_Ibip\tD"<<j<<"_Iam\tD"<<j<<"_Isoma\t";
    for(int j=1; j<=directions; j++)
        file<<"D"<<j<<"_V\t";
    file<<endl;
    */
    //cout<<"Middle cell: "<<(int)(AmCellLayer[0]->spikingChannel->units.size()/2)<<endl;
}

void AmacrineLayer::feedInput(const CImg<double>& image)
{
    CImg<double> temp(sx,sy,1,1,0.0);
    CImg<double> input(sx,sy,1,1,0.0);
    double xImg = 0.0;            // Direction X situated in a map centred on (xOffset,yOffset) of a amacrine cell
    double yImg = 0.0;            // Direction Y situated in a map centred on (xOffset,yOffset) of a amacrine cell
    double pixel_distance = 0.0;  // Distance of the reviewed pixel
    double pixel_angle = 0.0;     // Angle of the reviewed pixel
    double pixel_angle_gm = 0.0;  // Angle of the reviewed pixel for calculation of gamma multiplier
    double gammaMultiplier = 0.0; // Gamma multiplier for the Gaussian spatial filter in order of every dendrite's direction angle.
    double initial_ang = 0.0;     // Gives the initial angle of the input section for a specific dendrite
    double final_ang = 0.0;       // Gives the final angle of the input section for the same dendrite
    double size_ang = PI2/directions;
    double direction_ang = 0.0;   // For a given dendrite, variable Theta_zero to use in the spatial filter Gyam
    double xPos = 0.0;            // Variables that transform the x and y offset of a dendrite that are setted using the
    double yPos = 0.0;            // center of the axis on the middle of the image, different that on the real image
                                  // On an image, point (0,0) is on the superior left corner, so, it has to be changed.
    double area = M_PI*pow(radius,2);
    double areaPortion = area/directions;
    double iBip_norm_factor = 4*110.551890514714*directions*directions/area;
    double sxh = sx/2.0, syh = sy/2.0;
    //Enciso variables
    double h=0.0, sc1=0.0, sc2=0.0, H1=0.0, H2=0.0, dsc1=0.0, dsc2=0.0, gCl=0.0, v=0.0;
    int cell, direction;

    if(wTransient)
    {
      supplementaryTransient.feedInput( image );
      inputCurrent.feedInput( supplementaryTransient.read() );
    }
    else
      inputCurrent.feedInput( image );

    // Inputs for each somas
    for(int n=0; n<AmCellLayer[directions]->spikingChannel->units.size(); n++)
    {
        xPos = AmCellLayer[directions]->spikingChannel->units[n]->xOffset + sxh;
        yPos = syh - AmCellLayer[directions]->spikingChannel->units[n]->yOffset;
        v = AmCellLayer[directions]->spikingChannel->units[n]->get_V();
        for(int m=0; m<directions; m++)
        {
            input(xPos,yPos)+=delta*(AmCellLayer[m]->spikingChannel->units[n]->get_V()-v);
            file<<delta*(AmCellLayer[m]->spikingChannel->units[n]->get_V()-v)<<"\t";
        }
    }
    AmCellLayer[directions]->spikingChannel->feedInput(input,true);

    // For each direction, there's an input that must excite each dendrite in a separate way
      for(int n=0; n<directions ; ++n)
      {
          direction_ang = n*size_ang;

          if(n==0) initial_ang = PI2 - size_ang/2.0;
          else initial_ang = direction_ang - size_ang/2.0;

          final_ang = direction_ang + size_ang/2.0;

          input.fill(0.0);
          temp.fill(0.0);

          for(int m=0; m<AmCellLayer[n]->spikingChannel->units.size(); ++m)
          {
              temp.fill(0.0);
              xPos = AmCellLayer[n]->spikingChannel->units[m]->xOffset + sxh;
              yPos = syh - AmCellLayer[n]->spikingChannel->units[m]->yOffset;
              double inputs=0.0;

              for(double x = xPos - radius; x<= xPos+radius; x++)
              {
                  for(double y = yPos - radius; y<= yPos+radius; y++)
                  {
                      if(x < 0 || x >= sx || y < 0 || y>= sy) continue;
                      xImg = x - xPos;
                      yImg = - ( y - yPos );
                      pixel_distance = sqrt (pow(xImg,2)+pow(yImg,2));

                      pixel_angle = atan2 (yImg,xImg);
                      if (pixel_angle<0) pixel_angle = PI2 + pixel_angle;

                      pixel_angle_gm = atan2 (yImg,xImg);
                      if (pixel_angle_gm<0 && direction_ang!=0) pixel_angle_gm = PI2 + pixel_angle_gm;
                      if (pixel_angle_gm > (direction_ang+M_PI)) pixel_angle_gm -= PI2;
                      if (direction_ang>M_PI_2 && pixel_angle_gm<(direction_ang-M_PI)) pixel_angle_gm += PI2;

                      // Commented because the input receptors that receive from Bipolar cells cover all the dendrites area
//                      if(pixel_distance<=radius)
//                      {
                          // Because the first direction has angule 0 (n=0), the condition must be: higher than initial + PI OR lower than final
                          if( ( (n==0) && ( (pixel_angle>=initial_ang) || (pixel_angle<final_ang) ) )
                                  || ( (pixel_angle>=initial_ang) && (pixel_angle<final_ang) ) )
                          {
                              gammaMultiplier = (exp(pow((pixel_angle_gm-direction_ang),2)/(2.0*pow(gammaAm,2))) / (2.0*(M_PI)*pow(gammaAm,2)));
                              if(temp(x,y)==0)
                              {
                                  temp(x,y) = inputCurrent.read()(x,y) * gammaMultiplier * iBip_norm_factor;
                                  input(xPos,yPos) +=  temp(x,y);
                                  inputs++;
                              }
                          }
//                      }
                  }
              }
              input(xPos,yPos) *= areaPortion/inputs;
              file<<input(xPos,yPos)<<"\t";

              // Calculus about amacrine's interactions
              double distalC = 0.0;
              char outfile[200];
              snprintf(outfile,200,"../../../workspace/VR/SC2_%d.txt", n);
              //cout<<outfile<<endl;
              ofstream f, f2;
              f.open(outfile,ios::app);
              snprintf(outfile,200,"../../../workspace/VR/SC1_%d.txt", n);
              //cout<<outfile<<endl;
              f2.open(outfile,ios::app);
              h = fatherRetina->step;
              sc1 = AmCellLayer[n]->spikingChannel->units[m]->sc1;
              sc2 = AmCellLayer[n]->spikingChannel->units[m]->sc2;
              v = AmCellLayer[n]->spikingChannel->units[m]->get_V();
              H1 = AmCellLayer[n]->getSigmoidal(v, theta1, k1);
              dsc1 = alpha * ( 1.0 - sc1 ) * H1 - beta * sc1;
              sc1 = sc1 + h * dsc1;
              H2 = AmCellLayer[n]->getSigmoidal(sc1, theta2, k2);
              dsc2 = alpha * ( 1.0 - sc2 ) * H2 - beta * sc2;
              sc2 = sc2 + h * dsc2;
              f2<<sc1<<"\t";
              f<<sc2<<"\t";

              for(int j=1; j<AmCellLayer[n]->distalComp[m].size(); ++j)
              {
                  // Actual data
                  h = fatherRetina->step;
                  cell = AmCellLayer[n]->distalComp[m][j]->cellNumber;
                  direction = AmCellLayer[n]->distalComp[m][j]->direction;
                  sc1 = AmCellLayer[n]->distalComp[m][j]->sc1;
                  sc2 = AmCellLayer[n]->distalComp[m][j]->sc2;
                  //Calculus for sc1 for the new step
                  H1 = AmCellLayer[n]->getSigmoidal(AmCellLayer[direction]->spikingChannel->units[cell]->get_V(), theta1, k1);
                  dsc1 = alpha * ( 1.0 - sc1 ) * H1 - beta * sc1;
                  sc1 = sc1 + h * dsc1;

                  //Calculus for sc2 for the new step
                  H2 = AmCellLayer[n]->getSigmoidal(sc1, theta2, k2);
                  dsc2 = alpha * ( 1.0 - sc2 ) * H2 - beta * sc2;
                  sc2 = sc2 + h * dsc2;

                  AmCellLayer[n]->distalComp[m][j]->actualizeVariables(sc1,sc2);
                  distalC += sc2;
                  //cout<<j<<"\t"<<distalC<<endl;
                  f<<sc2<<"\t";
              }
              f<<endl;
              f.close();
              f2<<endl;
              f2.close();
              // Cl- channel conductance
              gCl = gclRest + (gclBound - gclRest)*distalC;

              // Adding the conductance input from amacrine cells to the dendrite
              AmCellLayer[n]->spikingChannel->units[m]->feedInput( ampInputCurrentAm*gCl,AmCellLayer[n]->amCondPort);
              file<<Ecl*gCl<<"\t";

              // Adding soma input
              v = AmCellLayer[n]->spikingChannel->units[m]->get_V();
              input(xPos,yPos) += /*areaPortion*/delta*(AmCellLayer[directions]->spikingChannel->units[m]->get_V()-v);
              file<</*areaPortion*/delta*(AmCellLayer[directions]->spikingChannel->units[m]->get_V()-v)<<"\t";
          }


          //Finally Ibip is given to the respective Amacrine Cells Spiking Channel
            AmCellLayer[n]->spikingChannel->feedInput(input,true);
            //AmCellLayer[n]->stock_input = input;
            temp = inputCurrent.read();
            AmCellLayer[n]->stock_input = temp;
            //cout<<temp.print()<<endl<<endl;
      }
}

void AmacrineLayer::tempStep()
{
    CImg<double> temp(sx,sy,1,1,0.0);
    double xImg = 0;            // Direction X situated in a map centred on (xOffset,yOffset) of a amacrine cell
    double yImg = 0;            // Direction Y situated in a map centred on (xOffset,yOffset) of a amacrine cell
    double pixel_distance = 0;  // Distance of the reviewed pixel
    double pixel_angle = 0;     // Angle of the reviewed pixel
    double initial_ang = 0;     // Gives the initial angle of the input section for a specific dendrite
    double final_ang = 0;       // Gives the final angle of the input section for the same dendrite
    double size_ang = 2*(M_PI)/directions;
    double direction_ang = 0;   // For a given dendrite, variable Theta_zero to use in the spatial filter Gyam
    double xPos=0;              // Variables that transform the x and y offset of a dendrite that are setted using the
    double yPos=0;              // center of the axis on the middle of the image, different that on the real image
                                // On an image, point (0,0) is on the superior left corner, so, it has to be changed.

    //for(int n=0; n<4 ; ++n)
      //  amacrineSignal[n].fill(0);

    stock_output.fill(0.0);

    if(wTransient)
        supplementaryTransient.tempStep();
    inputCurrent.tempStep();

    for(int n=0; n<directions ; ++n)
    {

        // The generated voltage Vam has to be placed at the right position on an image

        AmCellLayer[n]->stock_output.fill(0.0);
        direction_ang = n*size_ang;

        if(n==0) initial_ang = PI2 - size_ang/2.0;
        else initial_ang = direction_ang - size_ang/2.0;

        final_ang = direction_ang + size_ang/2.0;

        for(int m=0; m<AmCellLayer[n]->spikingChannel->units.size(); ++m)
        {

            //This "for" reviews every dendrite (for a given direction) and maintains its output section
            // for a specific direction of DSGC

            temp.fill(0.0);

            for(double x = xPos - radius; x<= xPos+radius; x++)
            {
                for(double y = yPos - radius; y<= yPos+radius; y++)
                {
                    if(x < 0 || x >= temp.width() || y < 0 || y>= temp.height()) continue;
                    xPos = AmCellLayer[n]->spikingChannel->units[m]->xOffset + ( temp.width()/2.0 );
                    yPos = (temp.height()/2.0 ) - AmCellLayer[n]->spikingChannel->units[m]->yOffset;
                    xImg = x - xPos;
                    yImg = - ( y - yPos );
                    pixel_angle = atan2 (yImg,xImg);
                    if (pixel_angle<0) pixel_angle = PI2 + pixel_angle;
                    pixel_distance = sqrt (pow(xImg,2)+pow(yImg,2));

                    if((pixel_distance>radiusOutputStarts) && (pixel_distance<=radius))
                        if( ( (n==0) && ( (pixel_angle>=initial_ang) || (pixel_angle<final_ang) ) )
                                || ( (pixel_angle>=initial_ang) && (pixel_angle<final_ang) ) )
                            temp(x,y) = AmCellLayer[n]->spikingChannel->units[m]->get_V();
                }
            }

            AmCellLayer[n]->stock_output += temp;
            stock_output += temp;
        }
    }

    for(int n=0; n<directions; n++)
    {
        for(int m=0; m<AmCellLayer[n]->spikingChannel->units.size(); m++)
//            if(m==(int)(AmCellLayer[n]->spikingChannel->units.size()/2))
                file<<AmCellLayer[n]->spikingChannel->units[m]->get_V()<<"\t";
    }
    for(int n=0; n<directions; n++)
    {
        for(int m=0; m<AmCellLayer[n]->spikingChannel->units.size(); m++)
//            if(m==(int)(AmCellLayer[n]->spikingChannel->units.size()/2))
                file<<AmCellLayer[n]->spikingChannel->units[m]->get_g_infini()<<"\t";
    }file<<endl;
}

AmacrineLayers::AmacrineLayers()
{
    outputSignal = 0;
    inputSignal = 0;
    amCondPort = 0;
}

AmacrineLayers::~AmacrineLayers()
{
    if(outputSignal)
        delete outputSignal;
    if(inputSignal)
        delete inputSignal;
    if(spikingChannel)
        delete spikingChannel;
}

double AmacrineLayers::getSigmoidal(double v, double theta, double k)
{
    return 1 / (1 + exp( -( v - theta ) / k ) );
}

void AmacrineLayers::addDistalComp(int baseCell, int dir, int cell)
{
    (this->distalComp[baseCell]).push_back(new distalCompartiment(dir,cell));
}

// ********************************************************
// ***************** DS Ganglion Layer ********************
// ********************************************************

DSGanglionLayer::DSGanglionLayer()
{
    sign=1.0;
    vThresh=0.0;
    valThresh=70.0;
    ampInputCurrent=1;
    amacrineAmpInputCurrent = 100;
    bipolarAmpInputCurrent = 100;
    tauTransient=0.04;
    wTransient=0.75;
    sigmaPool=0.0;
    layer = 0;
    IBip=0;
    nbCells=10000;
    inputSignal = 0;

    spikingChannel=0;
    gammaGang = 1;
}

void DSGanglionLayer::xmlize()
{
    add_param(sign,"sign");
    add_param(tauTransient,"transient-tau__sec");
    add_param(wTransient,"transient-relative-weight");
    add_param(vThresh,"linear-threshold");
    add_param(valThresh,"value-at-linear-threshold__Hz");
    add_param(amacrineAmpInputCurrent,"amacrine-amplification__Hz");
    add_param(bipolarAmpInputCurrent,"bipolar-amplification__Hz");
    add_param(sigmaPool,"sigma-pool__deg");
    add_param(nbCells,"number-of-cells");
    add_param(gammaGang,"gamma-direction-filter");

    add_child(spikingChannel,"");
    xmlParam::Register<1,ret_SpikingChannel>(spikingChannel,"spiking-channel");
    xmlParam::Register<1,ret_SpikingChannel_adap>(spikingChannel,"spiking-channel-adap");
}

void DSGanglionLayer::allocateValues(bool verbose, int chann)
{
    checkFather();

    int sx=fatherRetina->sizeX;
    int sy=fatherRetina->sizeY;
    double st=fatherRetina->step;
    double ppd=fatherRetina->pixels_per_degree;
    ret_LogPolarScheme *rs=fatherRetina->radScheme;
    double density=0.0, sxx=sx, syy=sy;
    layer += chann;
    if(verbose && !layer)
        cout<<"Allocating values for the selective-direction ganglion cells layer."<<endl;

    input = CImg<double>(sx,sy,1,1,0.0);

    inputSignal=new CImgReader;
    inputSignal->targets = &input;

    Iam = CImg<double>(sx,sy,1,1,0.0);

    inputAmSignal=new CImgReader;
    inputAmSignal->targets = &Iam;

    spikingChannel->set_step(st).set_ratiopix(ppd).set_radScheme(rs);
    spikingChannel->sampScheme->sizeX=&sxx;
    spikingChannel->sampScheme->sizeY=&syy;
    density = sqrt(nbCells/(sx * sy));
    spikingChannel->sampScheme->Density=&density;
    spikingChannel->sampScheme->set();
    spikingChannel->sampScheme->father = spikingChannel;
    spikingChannel->allocateValues();

    BaseRecFilter amacrineCuttingInGanglionar(st);
    amacrineCuttingInGanglionar.Exp(tauTransient).warningVerbose=false;
    amacrineCuttingInGanglionar*=wTransient;

    inputCurrent.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
    supplementaryTransient.set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);

    if(wTransient)
      supplementaryTransient-=amacrineCuttingInGanglionar;
    supplementaryTransient.allocateValues();

    inputCurrent.inputSynapse()
                .rectification(sign*vThresh,sign*ampInputCurrent,valThresh);
    inputCurrent.set_sigmaPool(sigmaPool);
    inputCurrent.allocateValues();

    if(fatherRetina->AmacrineCells)
        for(int n=0 ; n<fatherRetina->AmacrineCells->directions ; ++n){
            amSupplementaryTransient.push_back(new RetinaMapFilter);
            amInputCurrent.push_back(new RetinaMapFilter);

            amInputCurrent.back()->set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
            amSupplementaryTransient.back()->set_radScheme(rs).set_sizeX(sx).set_sizeY(sy).set_ratiopix(ppd).set_step(st);
            if(wTransient)
               *(amSupplementaryTransient.back())-=amacrineCuttingInGanglionar;
            amSupplementaryTransient.back()->allocateValues();

            amInputCurrent.back()->inputSynapse()
                        .rectification(sign*vThresh,sign*amacrineAmpInputCurrent,valThresh);
            amInputCurrent.back()->set_sigmaPool(sigmaPool);
            amInputCurrent.back()->allocateValues();
        }
}

void DSGanglionLayer::feedInput(const CImg<double>& vBip)
{
    if(wTransient){
        supplementaryTransient.feedInput(vBip);
        inputCurrent.feedInput(supplementaryTransient.read());
    }
    else inputCurrent.feedInput(vBip);

    if(fatherRetina->AmacrineCells)
    {
        if(wTransient)
            for(int n=0; n<(fatherRetina->AmacrineCells->AmCellLayer.size()-1) ; ++n){
                amSupplementaryTransient[n]->feedInput(fatherRetina->AmacrineCells->AmCellLayer[n]->outputSignal->read());
                amInputCurrent[n]->feedInput(amSupplementaryTransient[n]->read());
            }
        else
            for(int n=0; n<(fatherRetina->AmacrineCells->AmCellLayer.size()-1) ; ++n)
            amInputCurrent[n]->feedInput(fatherRetina->AmacrineCells->AmCellLayer[n]->outputSignal->read());
    }

}

void DSGanglionLayer::tempStep()
{
    double direction_ang = layer*(M_PI)/2;
    double amDirectionAng = 0.0;
    CImg<double> temp(fatherRetina->sizeX,fatherRetina->sizeY,1,1,0.0);

    input.fill(0.0);

    // Iam is passed through a spatial filter Gygang
    if(fatherRetina->AmacrineCells){
        double sizeAng = 2*(M_PI)/fatherRetina->AmacrineCells->directions;
        for(int n=0; n<fatherRetina->AmacrineCells->directions ; ++n)
        {
            amDirectionAng = n * sizeAng;
            amSupplementaryTransient[n]->tempStep();
            amInputCurrent[n]->tempStep();

            switch((int)layer) {

            case 0:
                    if(cos(amDirectionAng)<-0.001)
                    {
                        temp += cos(amDirectionAng) * amInputCurrent[n]->read();
                    }
                    break;
            case 1:
                    if(sin(amDirectionAng)<-0.001)
                    {
                        temp += sin(amDirectionAng) * amInputCurrent[n]->read();
                    }
                    break;
            case 2:
                    if(cos(amDirectionAng)>=-0.001)
                    {
                        temp += -cos(amDirectionAng) * amInputCurrent[n]->read();
                    }
                    break;
            case 3:
                    if(sin(amDirectionAng)>=-0.001)
                    {
                        temp += -sin(amDirectionAng) * amInputCurrent[n]->read();
                    }
                    break;
            default: cout<<"Non amacrine direction and DSGC layer could be related"<<endl;

            }

        }

        input += spikingChannel->feedInput(temp, gammaGang, direction_ang, (-1*amacrineAmpInputCurrent));
    }


    Iam = CImg<double> (input);

    if(wTransient)
      supplementaryTransient.tempStep();
    inputCurrent.tempStep();

    input += bipolarAmpInputCurrent * inputCurrent.read();

    // Now, Ibip is being charged
    spikingChannel->feedInput(input);
}
